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ABSTRACT 


An  Aerodynamic  Heating  and  Ablation  Computer  Program  (HEATAB)  Is  presented  to 
establish  a  means  by  which  heat  transfer  problems  may  be  solved  with  minimum 
effort.  This  program  computes  the  boundary  layer  conditions,  time-temperature 
distribution  in  a  body,  the  ablation  recession  rate,  and  weight  loss.  The 
program  was  written  in  Fortran  IV  for  the  CDC  6600  computer. 
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SECTION  I 
INTRODUCTION 

This  program  is  an  outgrowth  of  a  need  for  aerodynamic  heating  and  ablation 
analysis  of  hypersonic  reentry  vehicles.  The  present  design  requirements  for 
advanced  systems  of  high  performance  reentry  vehicles  have  resulted  in  a  severe 
thermal  environment  around  the  nose  and  conical  surfaces  of  the  vehicles.  The 
most  popular,  self-regulating  type  of  heat  protection  system,  the  ablating  heat 
shield,  will  produce  a  surface  erosion  and  nose  blunting.  Large  shape  changes, 
which  will  influence  the  aerodynamic  performance  of  the  vehicle  during  flight, 
may  be  produced.  This  vehicle  shape  change  for  an  ablation-type  system  is 
important  from  a  systems  standpoint  in  determining  the  actual  performance  of 
the  vehicle  during  flight.  Equally  important  from  a  systems  standpoint  is  the 
determination  of  boundary-layer  parameters  for  R-F  transmission  as  well  as 
temperature-time  profiles  in  the  vehicle  for  temperature-sensitive  instruments. 

It  is  our  hope  that  this  report  will  help  the  nonthermodynamicist  understand 
his  problem  and  solve  it  with  minimum  effort.  The  equations  of  heat  transfer 
are  general  in  that  they  apply  to  most  axially  symmetric  vehicles  of  any  compos¬ 
ite  skin  structure.  The  vehicle  skin  is  assumed  to  be  made  of  discrete  layers 
of  material  whose  properties  may  vary  from  layer  to  layer.  The  existing  com¬ 
puter  program  is  specific  in  that  the  equations  have  been  applied  to  a  sharp¬ 
nosed,  conical  vehicle.  With  a  thorough  understanding  of  thermodynamic  principles, 
as  applied  to  a  conical  vehicle,  and  with  a  basic  knowledge  of  computer  program— 
ming,  the  reader  should  be  able  to  adapt  this  program  to  any  vehicle. 
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SECTION  II 

BACKGROUND  AND  THEORY 

1.  Background 

Advances  in  hypersonic  atmospheric  flight  have  resulted  in  environments  cf 
extremely  high  temperatures.  Boundary  layer  temperatures  in  excess  of  10,000°? 
are  characteristic  of  vehicles  entering  the  atmosphere  at  hypersonic  speeds. 

As  a  result,  the  thermal  barrier  is  at  present  the  major  obstacle  to  the  safe 
reentry  of  vehicles.  The  term  "thermal  barrier"  refers  primarily  to  the 
problems  associated  with  the  dissipation  of  the  vehicle  energy  to  the  surrounding 
air.  For  sharp  bodies  at  hypersonic  speeds,  most  of  the  dissipation  occurs  in 
the  boundary  layer.  The  viscous  stresses  within  the  boundary  layer  exert 
shearing  work  on  the  fluid  and  raise  its  temperature  appreciably.  ‘This  work, 
called  aerodynamic  heating,  also  raises  the  surface  temperatures  of  bodies 
within  this  environment. 

In  hypersonic  flight,  it  becomes  obvious  that  capacitance  or  mass  heat  sink 
protection,  though  simple,  is  an  impossible  means  of  coping  with  the  extremely 
high  total  heat  input  associated  with  hypervelocity  reentry  vehicles.  Conse¬ 
quently,  other  means  of  cooling  or  protecting  for  operating  beyond  heat  sink 
limitations  must  be  used.  One  convenient  means  for  meeting  this  form  of 
environmental  protection  is  by  the  method  of  ablation  cooling. 

The  ablative  process  is  illustrated  in  figure  1.  The  ablation  process 
works  in  the  following  manner:  (1)  the  material  or  ablator  acts  as  a  heat  sink, 
(2)  when  the  critical  or  melting  temperature  of  the  ablator  is  reached,  a  thin 
layer  of  the  material  at  the  surface  will  begin  successively  to  melt,  vaporize, 
depolymerize,  or  decompose  chemically,  (3)  as  the  material  vaporizes,  the 
gaseous  products  of  decomposition  enter  into  the  boundary  layer.  Being  rela¬ 
tively  cool,  as  compared  to  the  boundary-layer  air,  the  injected  gas  forms  a 
thin  but  effective  film  that  greatly  reduces  the  heat  transfer  to  the  vehicle 
surface. 

Experiments  have  verified  that,  in  high-speed  flow,  the  magnitude  and 
direction  of  h  X  flow  at  the  surface  does  not  depend  on  the  difference  between 
the  wall  temperature  and  the  free-stream  temperature  as  in  low-speed  flow,  but 
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Figure  1.  Ablative  Process 

rather  on  the  difference  between  the  wall  temperature  and  the  adiabatic  wall 
temperature.  To  correlate  experimental  data,  it  is  convenient  to  define  the 
unit  surface  heat  rate  (reference  1)  for  high-speed  flow  as 

qc/A  -  he  (TAW  -  TW) 

where 

qc/A  *  heating  rate  BTU/sec  -  ft2 

he  *  convective  heat  transfer  coefficient  -BTU/ft2  -  sec  -  °R 
TAW  «  adiabatic,  wall  temperature  -  °R 
TW  *  vehicle  wall  temperature  -  CR 
A  =  local  wetted  area  -  ft2 

It  is  important  in  the  selection  of  a  suitable  material  for  thermal  protection 
by  the  ablative  process  that  the  products  of  decomposition  have  a  high  specific 
heat.  This  in  turn  produces  a  high  effective  mean  specific  heat  of  the  gas-air 
mixture  in  the  boundary  layer  and  a  high  Prandtl  number.  It  is  also  desirable 
that  the  ablator  have  a  low  thermal  conductivity.  This  will  decrease  heat 
conducted  to  the  interior  of  the  vehicle. 
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2.  Heat  Transfer  Coefficient 

In  conventional  analysis  of  heat  transfer  by  forced  convection,  the  nature 
cf  the  flow  can  be  described1 by  the  Reynolds  number,  which  is  a  dimensionless 
measure  of  the  ratio  of  inertial  to  viscous  forces.  In  high  speed  flow,  at 
least -two  additional  parameters  must  be  considered.  These  are  Mach  number  and 
Prandtl  number .  The  parameter  Mach  number  (AMACH)  describes  the  influence  of 
compressibility  on  heat  transfer  and  flow  phenomena  and  is  defined  as  the  ratio 
of  the  gas  or  flight  velocity  VEL  to  the  local  or  ambient  speed  of  sound  VS. 

The  Prandtl  number,  Pr,  is  defined  as  the  ratio  of  heat  storage  to  heat  conduc¬ 
tion  of  a  gas. 

The  heat  transfer  coefficient  (reference  2)  is  related  to  the  flow  properties 
through  the  relation 


Nu(Pr)*1/3  -  C(Re)a 


where 

Nu  -  Nusselt  number  «  he  x/K 

he  -  effective  heat  transfer  coefficient 

X  -  distance  from  the  nose  -  ft 

k  »  thermal  conductivity  of  air  -  BTU/ft  -  sec  -  °R 

Pr  «  Prandtl  number  *  C  u/K 

P 

Cp  -  specific  heat  of  air  -  BTU/lb  -  °R 

W  ■  viscosity  of  air  -  lb  -  sec/ft2 

Re  ■  Reynolds  nunber  *  pVx/p 

p  «  density  of  air  -  lb  -  sec2/ft4 

V  «  Boundary-layer  edge  velocity  -  ft /sec 

C  0,0296  for  turbulent  flow  (reference  1) 

0.354  for  laminar  flow 


a  *  0.8  for  turbulent  flow  (reference  1) 
0.5  for  laminar  flow 


Good  correlation  with  experimental  data  is  obtained  if  the  gas  properties  are 
evaluated  at  the  reference  temperature  and  the  velocity  is  taken  as  the  boundary- 
layer  edge  velocity.  For  the  case  of  mass  addition,  the  properties  of  the 
injected  gas  are  used  to  compute  the  Prandtl  number  and  Nusselt  number. 

3.  Unsteady  Heat  Conduction 

This  section  presents  a  means  of  applying  the  Schmidt  graphical  method  for 
solving  an  unsteady  heat  conduction  problem.  A  thorough  discussion  of  this 
method  can  be  found  in  any  general  text  on  heat  transfer  (references  1  and  3) 
and  will  not  be  attempted  here.  This  method  is  quite  flexible  and  difficult 
boundary  conditions  can  be  handled  easily. 

The  graphical  method  is  widely  used  in  industry  because  it  gives  a  running 
picture  of  the  changing  temperature  distribution;  its  details  can  be  delegated 
to  relatively  untrained  personnel,  and  mistakes,  if  they  occur,  are  quickly 
discovered.  However,  for  precise  computations,  especially  when  variations  of 
physical  properties  of  material  with  temperature  are  important,  as  in  the 
ablation  process,  a  numerical  method  should  be  used  instead  of  a  graphical 
method.  Numerical  methods  are  especially  convenient  when  a  computer  is  avail¬ 
able,  since  the  steps  involved  in  a  numerical  solution  can  be  programmed  without 
difficulty. 

The  numerical  method  for  solving  unsteady-state  conduction  problems  differs 
from  that  used  to  solve  steady-state  problems.  In  the  latter  cases,  the 
temperature  distribution  in  a  body  can  be  obtained  for  a  network  of  points  in  a 
solid  by  solving  a  system  of  residual  equations.  In  unsteady-state  systems, 
the  initial  temperature  distribution  is  known,  but  its  variation  with  time  must 
be  determined.  It  is,  therefore,  necessary  to  deduce  the  temperature  distribu¬ 
tion  at  some  future  time  from  a  given  distribution  at  an  earlier  time. 

To  illustrate  the  numerical  method,  it  is  first  necessary  to  transform  the 
Fourier  conduction  equation  (a  partial  differential  equation  that  is  second 
order  in  space  and  first  order  in  time)  for  the  unsteady  temperature  distribu¬ 
tion  in  a  heat-conducting  solid  into  a  finite  difference  form. 

AT  A2  T 

_T _ Ct  X 

At  ”  (Ax) 2 

the  subscript  denotes  the  differentiation  variable.  Letting  n  denote  position 
and  K  time,  A  T  can  be  written  as 
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\T  *  Tn,k+1  '  Tn,k 


In  a  similar  manner 


AxT  “  Tn+l,k  "Tn,k 


The  expression  A^2!  thus  becomes 


4x2t  •  W  -  2V +  W 


Substituting  these  expressions  into  the  Fourier  equation  gives 


T  .  -  T 

n,k+l  i 


:n,k  “  (Ax) 2  ^Tn+l,k  “  2Tn,k  +  Tn-l,k^ 


The  temperature  throughout  a  wall  or  slab  can  now  be  computed  for  any  later  time 
if  the  initial  distribution  is  known.  A  Schmidt  plot  demonstrates  the  tempera¬ 
ture  profile  versus  time  in  a  semi-infinite  slab  (figure  2). 

kTAW  - ORIGINAL  TEMPERATURE -#*0 

\  - TEMPERATURE  I  A#  LATER 

Y\  - TEMPERATURE  t  L I  LATER 

\\\  I  s  I  !  ! 


i , 

w. 


SK  ! 


hw 


sb?  iTt 

T* 


-AX — *f»— AX — — AX — AX — ■] 
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Figure  2.  Schmidt  Plot  in  an  Infinitely  Thick  Wall 


6 


A  constant  wall  temperature  is  assumed  in  the  graphical  example  but  a 
varying  wall  temperature  can  be  handled  with  equal  ease.  By  letting  the  wall 
temperature  vary  with  time,  the  distribution  at  each  point  within  the  body  can 
be  computed  for  each  time  increment. 

4.  Vehicle  Thermal  Model 


To  determine  the  time-temperature  profile  in  the  vehicle  skin,  the  Fourier 
conduction  equation  was  written  as  a  finite  difference  equation  and  solved 
numerically.  This  method  is  a  further  adaptation  of  the  Schmidt  graphical 
method  (reference  3).  The  wall  is  divided  into  a  number  of  lamina  with  known 
thickness,  specific  heat,  and  thermal  conductivity  and  with  a  known  initial 
temperature  distribution.  Figure  3  illustrates  the  model  with  its  associated 
nomenclature. 


AIK  FLOW 

- ► 


Figure  3,  Model  for  Wall  Temperature  Profile  and  Ablation 
Recession  Calculations 

The  difference  equation  to  compute  the  temperature  at  position  n  at  time 
k+1  is 


.k+1 


=  T 


h.k 


aAr 
(Ax)  2 


( 


n+l,k 


-  2T 


n.k 


+  Tn~l,k  ^ 
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In  this  fora,  the  properties  cannot  vary  from  layer  to  layer.  If  the  properties 
vary  from  layer  to  layer,  the  above  equation  must  be  written  in  the  form 


Tn,k+1  “  Tn,k  + 


-  V )  -  prV?  (v  - 


n+l,k 


It  is  useful  to  replace  the  thermal  diffusivity  by  its  defined  equivalent 


This  step  allows  one  to  easily  recognize  certain  groups  of  terms  as  the  heat 
conducted  through  each  lamina  and  permits  the  convective  heat  flux  to  be  used 
in  solving  for  the  wall  temperature.  The  equation  for  the  wall  temperature  is 


T  m  T  •f 

n-l,k+l  n-l,k  r 


^  [h‘(™ '  Vi,.)] 


f pcp  ** 


(w  -  v) 


The  backface  boundary  condition  can  be  specified  in  several  ways.  If  internal 
heating  (or  cooling)  is  present,  the  backface  temperature  may  be  held  constant 
or  allowed  to  vary  in  a  specified  manner.  A  conservative  assumption  is  that  no 
heat  flows  through  the  last  lamina.  This  assumption  will  cause  somewhat  more 
ablation  and  higher  temperatures  than  would  be  encountered  with  a  cooled  back- 


face  or  other  heat  sink  material. 


5.  Mass  Loss  and  Surface  Recession 


Two  mechanisms  are  involved  in  the  loss  of  ablation  material.  These 


mechanisms  are:  (1)  mass  loss  due  to  oxidation,  and  (2)  mass  loss  due  to 
sublimation,  The  oxidation  is  controlled  by  the  diffusion  of  oxygen  to  the 
reacting  carbon  and  the  sublimation  is  controlled  by  local  pressure  and 
temperature,  A  detailed  discussion  of  the  theory  of  ablation  material  and 
decomposition  of  reaction  is  not  the  intent  of  this  report,  so  only  the  results 
will  be  presented.  For  a  more  detailed  discussion,  see  references  4  and  5. 
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For  laminar  flow,  the  diffusion  mass  loss  rate  is  given  by 


•  _ qc  *  qHGR 

"D  *  K1t  +  K2t  I  hr  -  C_  TW 


L  L 


BL 


) 


and  for  turbulent  flow,  the  equation  is 


•  • 


“D 


qc  +  qHGR 


Kl„ 


+  K2T^ 


h  -  C  TW 


*  P 


BL 


) 


where  the  constants  K1  and  K2  can  be  interpreted  as  the  intercept  and  slope 
respectively  on  an  effective  heat  of  ablation  versus  enthalpy  difference 
(between  recovery  and  wall  conditions)  plot. 


q  =  heat  rate  -  BTU/ft2  -  sec 
c 

qFGR  “  heat  rate  ^or  hot  gas  radiation  “  BTU/ft2  -  sec 
h^_  =  recovery  enthalpy  in  boundary  layer  -  BTU/lb 
r 

PBL  =  sPec*f*c  beat  dn  boundary  layer  -  BTU/lb  -  °R 
TW  =  wall  temperature  -  °R 


qHCR  " 


GE  -  (TW)4  ALPG 


where 


o 

T 

GAS 

GE 

TW 

ALPG 

The  total 


-  Sf'phan-Boltzmann  constant,  4,8  x  10“ 13 
=  gas  temperature  -  °R 
=  gas  emissivity  at  temperature  T^^, 

=  wall  temperature  -  °R 
=  absorptivity  of  gas  at  temperature  TW 
mass  loss  rate  is  given  by 


BTU/sec-ft2  - 


°R4 
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V*D 

where 

PBL  *  bo’mdary-layer  edge  pressure  -  lbs /ft2  (reference  6) 
The  surface  recession  rate  is  calculated  by 


1  +  2.64  x  104(PBL)“°*67  exp 


(- 


11.05  x  104 
TW 


) 


S  -  M^p 

where  p  is  the  ablation  material  virgin  density.  To  determine  the  total 
recession,  S  is  summed  after  each  calculation  for  a  new  trajectory  point. 

6.  Calculation  of  Total  Weight  Loss  Rate 

After  the  local  instantaneous  and  tc*al  recession  rates  have  been  computed, 
the  local  weight  loss  rate  is  calculated  by  determining  the  volume  reduction  of 
the  body  and  multiplying  by  the  material  density.  It  is  assumed  in  this  calcu¬ 
lation  that  the  recession  is  uniform  over  the  region  of  interest.  Total  weight 
loss  rate  is  obtained  by  summing  the  local  rates. 
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SECTION  III 
PROGRAM  DESCRIPTION 

1.  Introduction 

This  section  discusses  ana  asents  the  basic  program  (HEATAB)  with  stagna¬ 
tion  point  modification  (STAGG)  for  a  sharp-nose,  conical  vehicle.  The  program 
was  written  in  FORTRAN  IV  for  the  CDC  6600  computer.  Only  a  general  flow  chart 
is  given.  The  flow  within  individual  blocks  is  easily  followed  from  the 
program  listing  (Appendixes  I  and  II) . 

2.  HEATAB  Description 

The  basic  sequence  of  calculations  is  shown  on  the  flow  chart  of  HEATAB, 
figure  4.  After  entry  into  the  program,  the  first  steps  are  the  reading  of 
input  data  and  the  defining  of  constants.  Descriptive  initial  conditions  are 
printed  to  identify  the  case  being  run. 

Time  is  the  independent  variable  and  its  computation  marks  the  entry  to  the 
computational  loop  of  the  program.  The  flight  conditions,  atmospheric  properties, 
and  layer  index  are  now  determined.  A  number  of  methods  can  be  used  to  input 
the  flight  trajectory  such  as  subroutine,  punched  card,  or  magnetic  tape.  The 
atmospheric  properties  are  computed  by  subroutines.  K  is  a  subscript  that 
indicates  the  layer  for  which  computation  is  being  performed.  The  initial 
value  for  K  is  1.  IND(K)  is  used  to  indicate  to  the  computer  to  use  virgin  or 
char  properties.  If  T(K)  is  less  than  the  char  temperature,  IND(K)  will  route 
the  computer  to  virgin  material  physical  properties.  Similarly,  if  T(K)  is 
greater  than  the  char  temperature,  then  char  properties  will  be  computed. 

The  next  sequence  of  calculations  determines  the  recovery  temperature, 
boundary-layer  edge  conditions,  and  Reynolds  number.  Following  these  calcula¬ 
tions,  it  is  determined  if  ablation  is  occurring  or  has  occurred.  Based  on  this 
determination,  the  correct  Prandtl  number  is  computed.  Based  on  the  outcome  of 
the  decision  regarding  the  existence  of  turbulent  flow,  the  appropriate  recovery 
factor,  heat  transfer  coefficient,  dynamic  temperature,  adiabatic  wall  tempera¬ 
ture,  and  gas  temperature  are  calculated.  The  resulting  output  of  this  sequence 
is  the  heat  flux  to  the  body. 
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Computing  the  wall  temperature  distribution  is  the  function  of  the  next 
sequence  of  computations.  Again  the  decision  regarding  charring  and  ablating  is 
made  and  the  appropriate  physical  properties  are  selected.  The  temperature 
TW(1)  is  computed  as  a  function  of  the  adiabatic  wall  temperature  and  TW(2) . 

The  temperatures  below  the  surface  are  next  computed  as  functions  of  the  old 
adjacent  temperatures. 

If  the  wall  is  ablating,  mass  is  being  lost  from  the  vehicle  and  the  wall 
will  recede.  The  mass  loss  rate  is  computed  in  this  sequence,  along  with  the 
recession  rate  and  new  layer  thickness.  The  new  thickness  is  used  the  next 
time  the  temperature  calculations  are  made.  The  printing  of  output  data  marks 
the  completion  of  one  cycle  of  computation.  As  a  new  cycle  is  started,  a 
decision  is  made  based  on  the  vehicle's  altitude.  If  the  altitude  for  the  new 
time  is  positive,  computation  continues;  if  the  altitude  is  negative,  a  normal 
exit  is  made  and  computation  ends. 

3.  STAGG  Description 

This  section  contains  a  description  of  the  STAGG  version  of  the  aerodynamic 
heating  and  ablation  program.  STAGG  is  a  modification  of  the  basic  program 
that  can  be  used  to  calculate  heat  transfer  rates  on  the  stagnation  point  of 
sharp-nosed  vehicles.  Heat  transfer  rates  and  temperature  profiles  are  calcu¬ 
lated  by  the  existing  HEATAB  equations,  using  local  properties,  but  the  constant 
in  the  heat  transfer  coefficient  and  Reynolds  number  for  transition  from  laminar 
to  turbulent  flow  has  been  redefined.  The  constants  as  defined  by  Lees'  theory 
(reference  7)  and  Blausius'  skin  friction  coefficient  for  point  cone  heat  flux 
are 

C  »  0.778  for  laminar  flow 
0.0348  for  turbulent  flow 

Another  difference  between  HEATAB  and  STAGG  is  the  addition  of  nose  blunting 
equations. 

In  the  determination  of  nose  blunting,  it  is  assumed  that  the  initial  nose 
shape  of  particular  interest,  for  which  STAGG  was  written,  was  a  sphere-cone 
and,  after  ablation,  the  final  shape  was  a  sphere-cone.  Experimental  evidence 
has  shown  that  the  final  eroded  shape  could  be  approximated  by  a  sphere-cono 
having  a  small  deviation  from  a  sphere.  The  reference  sphere-cone  configur  .ion 
after  erosion  is  obtained  by  placing  a  spherical  surface  tangent  to  a  cone 
parallel  to  the  original  surface  but  displaced  by  the  erosion  on  the  cone  and 
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and  located  axially  at  a  position  on  the  vehicle  center  line  in  accordance  with 
the  erosion  at  the  stagnation  point.  (See  figure  5.) 

The  experimental  evidence  of  nose-shape  change  of  reentry  vehicles,  as 
previously  noted,  can  be  represented  by  a  sphere-cone  shape  which  deviates 
from  a  spherical  shape  by  an  amount  less  than  10  percent  of  the  final  nose 
radius.  The  final  eroded  nose  radius  at  zero  angle  of  attack,  therefore,  can 
be  obtained  in  terms  of  the  stagnation  point  erosion  and  cone  erosion  as  shown 
in  equation 


R  -  R  +  — : — ~ — 
o  (1-Pin9c) 


sin0 


f:  (k)  "■/:($  ‘ 


The  integral  represents  the  erosion  at  a  given  time  during  a  flight  for  the 
locations  at  the  stagnation  point  (S)  and  a  point  on  the  cone  (C)  at  a  location 
of  wetted  length  of  about  five  times  the  nose  radius. 

It  is  now  possible  to  evaluate  the  final  nose  radius  for  a  given  application. 

Let 


i 


*  .R  =  R  +  — ~ —  T  sin6  TDOT  -  CDOT  1 

o  (l-sin9c)  [  c  J 


where 

Rq  =  original  nose  radius  -  ft 
R  =  instantaneous  nose  radius  --  ft 
9c  =  cor.e  half-angle  -  radians 


% 
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TDOT  =  instantaneous  stagna  Ion  point  recession  rate  -  ft/sec 
CDOT  =  instantaneous  cone  recession  rate  -  ft/sec 
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SECTION  IV 

v 

BOUNDARY  LAYER  TRANSITION 

The  problem  of  transition  from  laminar  to  turbulent  boundary  layer  flow  has 
always  been  a  troublesome  one  for  the  aeronautical  engineer.  The  standard 
approach  has  usually  been  to  design  conservatively,  that  is,  for  turbulent  flow. 

In  general,  the  transition  Reynolds  number  has  been  found  to  depend  upon  the 
local  Mach  number,  which  has  been  observed  in  the  analysis  of  boundary  layer 
transition  on  a  series  of  flight  vehicles.  Transition  on  blunt  spherical  nose 
vehicles  with  low  edge  of  boundary  layer  Mach  numbers  appears  to  occur  at 
Reynolds  number  of  1.5  x  105,  while  on  sharp  bodies  where  the  edge  of  the 

boundary  layer  Mach  numbers  approach  10,  Reynolds  numbers  as  high  as  2.0  x  106 

\ 

have  been  observed  prior  to  transition.  Thus,  for  untested  vehicle  configurations, 
it  is  necessary  to  investigate  the  relationship  of  the  Mach  number  and  the 
Reynolds  number  in  both  the  high  and  low  Mach  number  regions.  It  should  be 
noted  that  flight  data  on  sharp  bodies  indicates  that,  generally,  transition 
does  not  occur  instantaneously  over  the  entire  vehicle  so  that  it  may  travel 
several  thousand  feet  between  the  onset  of  transition  and  the  establishment 
of  a  fully  turbulent  boundary  layer.  These  criteria  have  been  summarized  by 
Hm.i;  in  reference  8. 

The  approach  of  this  report  is  to  apply  the  sharp  and  blunt  body  transition 
Reybolds  numbers  simultaneously.  Therefore,  the  sharp  body  criterion  is 
described  by  the  expression 

Re_  -  2.0  x  106 
tran 

where  the  Reynolds  number  is  that  from  a  wetted  length  of  one  foot  to  the 
extreme  aft  end  of  a  conical  vehicle. 

The  blunt  body  criterion  is  applied  when  either  of  the  following  two 
Reynolds  numbers  are  reached: 

Retran  =  1.5  x  105  on  the  spherical  nose 

Retran  =  x  on  conical  portion  aft  of  the  spherical  nose  for 
a  wetted  length  of  five  times  the  nose  radius 

The  above  criteria  has  been  found  to  correlate  with  experimental  data. 
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APPENDIX  I 

HEATAB  Program  Listing 
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LIST  OF  SUBROUTINES* 


Subroutine 

Computes 

Versus 

ALTUDE 

Alt,  Altitude  of  Vehicle 

Time 

CCPP 

CC,  Char  Specific  Heat  of  Lamina 

K  of  Ablator 

TW(K) 

CTCPC 

CT,  Char  Thermal  Conductivity 
of  Lamina  K  of  Ablator 

TW(K) 

DENS 

Rho,  Density  of  Air 

Alt 

GASE 

GE,  Ablative  Gas  Emissivity 

TG 

GASTC 

GTC,  Ablative  Gas  Thermal  Conductivity 

GCPB 

GCP,  Ablative  Gas  Specific  Heat 

P 

PRESS 

PE,  Atmospheric  Pressure 

Alt 

QVST 

Q,  Dynamic  Pressure 

Time 

SPCh 

SPC,  Cone  Surface  Pressure  Coefficient 

AMACH 

TCAT 

TCA,  Thermal  Conductivity  of  Air 

TR 

TEMPA 

Temp,  Ambient  Temperature 

Alt 

VCPPC 

VC,  Virgin  Specific  Heat  of  Lamina  K 
of  Ablator 

TW(K) 

VELOC 

Vel,  Velocity  of  Vehicle 

Time 

VI  SCO 

Vis,  Viscosity  of  Ambient  Air 

Alt 

VI SCOT 

VISG,  Ablative  Gas  Viscosity 

TR 

VSD 

VS,  Velocity  of  Sound 

Alt 

VTCPC 

VT,  Virgin  Thermal  Conductivity  of 

Lamina  K  of  Ablator 

TW(K) 

*This  listing  of  subroutines  common  for  HEATAB  and  STAGG, 
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ACP(K) 

AJ 

AKL1 

AKL2 

AKT1 

AKT2 

ALPG 

ALT 

AMACH 

ANR 

ARHO(K) 

ATC(K) 

CPA 

DE(K) 

DEL(K) 

DELT 

DMD 

DTR 

E 

G 

GCP 

GTC 

GVIS 

HE 

IND(K)  -  0 

IND(K)  «  1 

JIND  *  0 

JIND  »  1 

NAL 

NPl 

NP2 

P 

PE 


LIST  OF  SYMBOLS  FOR  HEATAB 

Specific  heat  of  ablator  (char  or  virgin),  BTU/lb  -  °R 
778.184,  Joule’s  constant 

Ablator  constants  for  mass  loss.  Can  be  interpreted  as  the 
intercept  and  slope  respectively  on  an  effective  heat  of 
ablation  versus  enthalpy  difference  plot  for  laminar  flow. 

Same  as  above,  except  for  turbulent  flow 

Absorptivity  of  ablative  gas 
Altitude  of  vehicle,  ft 
Mach  number  of  vehicle 
Computed  Reynolds  number 

Density  of  ablator  at  lamina  K  (char  or  virgin),  lbs/ft3 

Thermal  conductivity  of  ablator  at  lamina  K  (char  or  virgin) , 
BTU/ ft-sec-°R 

Specific  heat  of  air,  BTU/lb-°R 
Instantaneous  lamina  thickness,  ft 
Original  lamina  thickness,  ft 
Time  step,  sec 

Diffusion  controlled  oxidation  mass  loss 

Dynamic  temperature  rise,  °R 

2.7183 

Gravitational  acceleration,  ft/sec2 

Ablation  gas  specific  heat,  BTU/lb-°R 

Ablation  gas  thermal  conductivity,  BTU/ft-°R-sec 

Boundary  layer  air  viscosity,  lb-sec/ft2 

Convective  heat  transfer  coefficient ,  BTU/ft2-sec-  °R 

Lamina  K  of  ablator  is  virgin 

Lamina  K  of  ablator  is  char 

Outermost  lamina  not  ablating 

Outermost  lamina  is  ablating 

Number  of  ablative  laminas 

Number  which  refers  to  bond  lamina 

Number  which  refers  to  aluminum  skin 

Boundary  layer  pressure,  Atmospheric 

Ambient  pressure,  psi 
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PBL 

PR 

PXE 

Q 

QDOTC 

QDOTR 

RE 

RF 

RHO 

RHOC 

RHOV 

S 

SDOT 

SIGMA 

SPC 

TAW 

TCA 

TEM 

TEMP 

TG 

TIME 

TM 

TMDOT 

TMELT 

TR 

TW(K) 

VEL 

VIS 

VISG 

VS 

X 

WT 


Boundary  layer  edge  pressure,  psi 

Prandtl  ntanber 

Dummy  variable 

Dynamic  pressure,  psf 

Convective  heat  rate,  BTU/ft2-sec 

Hot  gas  radiation  heat  rate,  BTU/ft2-sec 

Transition  Reynolds  number 

Recovery  factor 

Ambient  density  of  air,  lb-sec2/ft3 

Char  ablator  density,  lb/ft3 

Virgin  ablator  density,  lb/ft3 

Instantaneous  ablator  wall  recession  rate 

Total  ablator  wall  recession,  ft 

Stephan-Boltzmann  constant 

Cone  surface  pressure  coefficient 

Adiabatic  vail  temperature,  °R 

Thermal  conductivity  of  air,  BTU/ft-°R-sec 

Initial  temperature  of  laminas,  °R 

Ambient  temperature,  °R 

Ablation  gas  temperature,  °R 

Elapsed  time,  sec 

Total  mass  from  ablation 

Instantaneous  mass  loss  rate 

Temperature  at  which  ablation  starts,  °R 

Reference  temperature,  °R 

Temperature  of  lamina  K,  °R 

Tangential  velocity  of  vehicle,  ft/sec 

Viscosity  of  ambient  air,  lb-sec/ft2 

Ablative  gas  viscosity,  lb-sec/ft2 

Velocity  of  sound,  ft/sec 

Wetted  length  from  nose  to  station  on  cone,  ft 
Instantaneous  weight  loss  rate 
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PROGRAM  HEATARI INPUT .OUTPUT) 

0  AERODYNAMIC  HEATING  ON  CONE 

DIMENSION  TWI50) . INDI50) , ATCI50) , ACPI  50) .ARH0I50) ♦ DEL (50)  .DFI50) 
READ  1 »  A  J  »  AKL 1  , AKL2  «  AKT 1 , AKT2 . G • GV I S * ? I GMA , NAL . ND 1 , NP2 , E , TMELT , 
1ARH0INP1  ) «  ARHO ( ND2 ) «  ATC I NP 1  ) .  ATC I NP2 ) . ACP(NP1 ) *ArP(  NP? )  «  DEL ( NP 1  )  , 
2DELINP?) ,DE(NPl ) ,Df (NP2) , DELT , RE , VEL , TEMP , CPA , RHOV , 

3RH0C  « TEM  «  TOL  *  X 

I  FORMAT!  AS  NECESSARY  ) 

DO  P  N= 1 «NP2 

TWIN) =TFM 
?  I ND ( N ) =n 

DO  3  N=1  . NAL 
DEL  IN) =TOL 

3  DE(N)=TOL 
SDOT  =  0  • 

JIND=0 
TIME=-DFLT 
TAW=TEMP 
PRINT  4 

4  FORMAT!  HEADINGS  AS  NECESSARY  ) 

100  T I Mt=T I ME+DELT 

DO  6  N= 1  « NAL 
AN=N 

IF  ( AN*DEL I  1  )  -  SDOT)  6,8.8 

6  CONTINUE 
PRINT  7 

7  FORMAT  I  1  OH  OVFRHEAT / ) 

GO  TO  200 

8  K  =  N 

TR=  TEMP+ . 58* I T W I K ) -TEMP )  +  . 1 Q* ( TAW-TEMP ) 

CALL  ALTUDE I ALTf TIME  ) 

CALL  VSO(VS«ALT) 

CALL  VELOC I VEL , T I  ME ) 

AMACH=VEL/VS 

CALL  TFMPA I  TEMP, ALT ) 

CALL  RHOF IRHO, ALT  ) 

CALL  V  I  SCO  I  V  I S  » ALT ) 

CALL  TCAT(TCA,TR) 

CALL  V  I  SCOT  I  V  I SG , TR ) 

CALL  SPCM ( SPC , AMACH ) 

CALL  OVS f I Q , T I  ME ) 

CALL  PRESSIPE, ALT) 

PPL  =  SPC  *0+PF 
P  =  P PL/21  1 6.2 
I F ( J I ND .EG • 0 ) GO  TO  Q 
CALL  GCPRIGCP.P) 

CALL  GASTC ( GTC ,GCP ) 

PR=GCP*GV I S/GTC 
00  TO  10 

o  pQ=CPAtty 1  SG/TCA 

10  ANP=VEL*X*PHO/V I c 

I^I ANP.GT.RF )G0  TO  12 

I I  PF=PQ**,5 

Hr  =  . 575*TCA/X*PR**.  33*  A  NR**.-. 

GO  TO  13 

1 2  or  =  pp**,3  3 
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HF  = .02Q6*TCA/X*PR**. 33* A  NR** . 8 
17  DTR=VEL**2./< 2.*G*AJ*CPA ) 

T  Aw=TFMP+PF  *DTP 

TG=TE^P+OTR 
CALL  GAFF ( GE ,  TG ) 

ALPG= • 022* ( TG/TMELT ) #*  «  45 
ODCTC=H£*<TAW-TW(K) )*DELT 

14  I F (  I  NO (<) •  EO  • 1 )GO  TO  17 
IF(TW(K) .GF.TMELTIGO  TO  17 

15  TWK=TW(K) 
call  VTCPC( VT,TWK) 

CALL  VCPOC( VC,TWK) 

ATC ( K ) = VT 
ACP ( < ) = VC 
ARHO<K) =RHOV 
IF ( k • FQ  # 1 ) GO  TO  18 

TU'(K)=TW(K)+DELT*<ATC<K-1  >/(ARHO(K-l  )*ACP(K-1  )*DE  (K-l  )**?)* 

1  (TW1K-1  )-TW(K)  )-ATC<K)/< ARHO < K ) *ACP < K > *DE  ( K ) **2 ) * ( TW ( K ) -TW ( K+ 1  )  ) 
1F(TW(K) .LT.TMELT5GO  TO  31 

16  TW<K)=TMFLT-1 • 

GO  TO  "M 

17  TWK=TW<K) 

CALL  C.TCPC(CT,TWK) 

CALL  CCPP ( CC  «  TWK ) 

ACP(K) =CC 
ATC(K) =CT 
ARHO { K ) -RHOC 
I  NO (<)  =  1 

18  TW(K ) =TW(K)+DFLT/< ACP ( K ) *ARHO < K ) *DE  ( K )  ) * ( HE* ( T AW-TW ( <  )  ) - 
1 ATC ( K ) /PE  { K ) * ( TW ( K ) — TW ( K+ 1 ) ) ) 

IP  I F ( T W ( K ) • GF • TMELT ) GO  TO  81 

20  jind  =0  r 

GO  TO  31 

21  TW ( K ) =TMELT 
I  NO ( K ) s 1 
J I N0=  1 

CALL  GCPP ( GCP « P 1 

QPOTR=F  I  GMA*  (  GE*TG**4-ALPG*TMELT**4  )  *DELT 
IF(0D0TQ)22,22,23 

22  OnOTPaO. 


\" 


V  }  \  i 

V  *  *  v V •  \ 

’s  .  I 


27 

24 

25 

26 


2P 

30 


1F< ANR.GT.PFJGO  TO  25 

DMD= ( QPOTC  +  ODOTR  > / ( AKL 1 +AKL2*  <  GCP*T AW-CPA*TMFLT )  ) 
GO  TO  26 

DMD= { ODOTC+QDOTR ) / ( AKT1 +AKT2* ( GCP*T AW-CPA*TMFLT ) ) 
PXE  = 1  1  .05F+04/TMELT 

TMDOT  =  DMD*{ 1 •  +  2 . 64E+0P/ < PBL** . 67*E#*PXE ) ) 
tm=tm+tmdot 

FOOT  =  TM/RHOV 
F=TMPOT/RHOV 

WT  =  3 • 14 16*2»*R**2*S*PH0V 
DO  2P  Mr  1 ,NAL 
5N  =  N 

!F  <BN*nFL<l>  -  FOOT)  2P, 30,30 
CONTINUE 
dp  j  MT  -r 
GO,  TO  200 

OF ( N ) =nFL ( M ) - < FOOT- ( N- 1 ) *DEL ( 1) ) 

25 


i 

f 

\ 


I NO ( N ) a  1 

IFfN.EO.llGO  TO  31 
DE ( N— 1 1*0.0 
31  K*K+1 

IF(K.UE.NAL1G0  TO  14 

TW(K)*Tw(K1+0ELT*(ATCCK-1  )/( ARH0(K-1  )*ACP(K>1  >*DE(K-1  )**2>*(TW<s<. 
l-ll-TWOO  >-ATC(K)/tARH0<t<)*ACP<K)*DFLCK:>**2)*(TW<K)-TW(K+l  )  )  ) 
K*K+l 

TW(K)*TW<K)+0ELT*(ATC(K-1  )/<ARHO<K-l  >*ACP(K-1  )*DEL(K-1  )  **2  )  * 

1 (TW1K-1 )-TW(K) 1 1 

PRINT  5.TIME.ALT.VEL.ANR,TAW.TEMP,HE,QD0TC«WT.TM.TG.S00T. 

1 (TWIN) «N=t  «NP2) 

5  FORMAT (  AS  NECESSARY  > 

IF ( ALT 1200.200.100 
200  CONTINUE 
ENO 
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CDOT 

LIST  OF  SYMBOLS  FOR  STAGG 

Instantaneous  cone  recession  rate,  ft/sec 

PHI 

See  figure  5 

R 

Radius  of  spherical  nose,  ft 

TCAE 

Time  cone  ablation  ends 

TCAS 

Time  cone  ablation  starts 

THETC 

Cone  half-angle,  rad 
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PPOGP AM  <;t  AGG  (  I NPUT  ,  OUTPUT  ) 

C  STAGNATION  HEATING 

DIMENSION  TVM50) , INDI50) ,ATC<50) « ACP < 50 ) . ARHO ( 50 ) » DEL (50) «DE (50) 
READ  1  «  AJ«AKL1  «  AKL2  ,  AKT  1  «  AKT2  *  G,  GV  I  S  *  $.  I  GMA  ,NAL  »NP1  ,  NP2  ,  E  ,  TMELT  , 
1ARH0(NP1 ) .  ARHO ( NP2 ) ,ATC<NPl ) ,ATC(NP2 ) ,ACP(NP1 ) ,ACP(NP2) ,DEL(NP1  )  * 
2DEL  <  NP2 ) ,DE(NP1 ) ,DE<NP2) ,DELT , P . THETC , PH  I ,RE « VEL , TEMP » CPA  ,  RHOV « 
3PH0C  ,  TO AS  , TEM , TOL , TC AF 

1  FORMAT <  AS  NECESSARY  ) 

DO  2  N=I « NP2 
TW(N)=TEM 

2  IND(N) =n 

DO  3  N=1 »NAL 
DEL(N) =TOL 

3  DE<N)*TOL 
CDOTaO.O 
SD0T=0. 

JIND=0 

T!ME=-OFLT 

X=.0t745*R*PHl 

TAw=TEMP 

op  I  NT  4 

4  FORMAT (  HEADINGS  AS  NECESSARY  > 

100  T I ME=T I MF+DELT 

00  6  N=1 « NAL 
AN=N 

IF  <AN*DEL(1)  -  SDOT )  6.8,8 

6  CONTINUE 
PRINT  7 

7  FORMAT (1  OH  OVERHEAT/) 

GO  TO  200 

R  K*N 

TR  =  TEMP+.58*<  Tw(K)-TEMP)  +  . 19* (TAW-TEMP) 

CALL  ALTUDE< ALT. TIME) 

CALL  VSD(VS.ALT) 

CALL  VFLOC( VEL, TIME) 

AMACHrVFL/VS 

CALL  TFMPA( TEMP, ALT) 

CALL  RHOF(RHO,ALT> 

CALL  V  I  SCO ( VIS, ALT ) 

CALL  TCAT(TCA,TR) 

CALL  V  I  SCOT ( V  I ^G , TR ) 

CALL  SPCM(SPC,AMACH) 
call  0VST(0,TIME) 

CALL  PPP^S ( PE , ALT ) 

PBL=SPC*Q+PE 
P=P8L/2 1 16.2 
I F ( J I ND  «  FQ , 0 ) GO  TO  9 
CALL  GCOR(GCP,P) 

CALL  GAOTC(GTC,GCP) 

PP=GCP*GV?  S/GTC 
GO  TO  10 

Q  PR=CPA*V I SG/TCA 

10  ANP=VEL*Y*PHO/VI  «=■ 

I F ( ANP  » GT ,PE ) C 0  TO  1? 

1 1  RF=PQ**,^ 

HF= ,77R*TCA/X*PP**. 37* ANP**, 5 
GO  TO  1 3 


DF:PO*#,?T 

HE=.0296*TCA/X*PR**. 33*ANR**. 8 
DTR  =  VEL**2./<  2.*G*AJ*CPA ) 
TAW=TEMP+RF*DTR 
T  G=  TEMP+PTR 
CALL  GAGtr(GF«TG) 

ALPG= • 022* ( TG/TMELT ) ** .45 
GD0TC=HP*( TAW-TW(K) )*DELT 
1 E (  I ND  (K).FO.l) GO  TO  17 
IF(TW<K)  .GE. TMELT) GO  TO  17 
TWK=TW(K> 


CALL  VTCPC( VT.TWK) 

CALL  VCPPC.  (  VC  «  TWK  ) 

ATC (K )  =V’T 

ACP(K) =VC 

ARHO ( K ) =PHOV 

IF ( K . EO « 1 ) GO  TO  18 

TW( K ) =TW ( K ) +DELT* ( ATC (K  -  1  )/<ARHO(K-l  )*ACP(K-1  )*DE  (K-l )**2>* 

1  <TW(K-1 )-TW(K) ) — ATC ( K ) / ( ARHO ( K ) *ACP ( K ) *DE  ( K ) **2 ) * ( TW ( K ) -TW ( K+ 1  )  )  ) 
I F ( TW ( K ) »LE .TMELT ) GO  TO  31 
TW(K)=TM<^LT-1 . 


GO  TO  31 
TWK=TW  (K  ) 

CALL  CTCPC(CT.TWK) 

CALL  CCDP( CC, TWK ) 

ACP ( K ) =CC 
ATC(K)=CT 
ARHO(K)=RHOC 
I ND  <  K )  =  1 

TW(K)=TW(K)+DELT/< ACP(K)*ARHO(K)*DE  <  K  >  ) * ( HE* ( T  AW— TW ( K )  )- 
1 ATC ( K ) /DE  (K)*(TW(K>-TW<K+1  )  )  ) 

I F ( TW ( K ) »GE. TMELT )G0  TO  21 
JIND  =  0 
GO  TO  71 
TW ( K ) =  TMFLT 
I ND (  K  )  =  1 
J  I  ND=  1 

CALL  GC.PP  (  GCP  •  P  ) 

ODOTR=c I GMA* ( GE*TG**4-ALPG*TMELT**4 ) *DELT 
I  F  (  QDOTP  )  22  » 2? 


ODOTR=0. 


IF(ANR.GT.R,  25 

DMD= ( ODOTC  +  ODOTR ) / <  AKL 1 +AKL2* ( GCP*T AW~CPA*TMELT  )  ) 
GO  TO  26 

DMD= ( ODOTC  +  OOOTR ) / ( AKT 1 + AKT2* ( GCP*TAW-CPA*TMELT  )  ) 
PXF  = 1  1 .P5F  +  04/TMFLT 

TMDOT  =  DMD*(1.  +  2  * 64F  +  09/ ( PBL** .67*F**PXE ) > 

TM=TM+TMDOT 

FOOT  =  TW/RHOV 

«=TMD0T/PHOV 

WT=3. 1416*?.*R**2*^*RH0V 

I  F  (  T  !  MF »LT  »  TC  A  c  .OR.  T I  ME . GT . TC AE ) GO  TO  28 
RF AP  27, COOT 
PQRMAT ( F 1 1 .5 ) 

R=R+ i ./(  1  ,-SI N< THETC )  )*( SIN<  THETC )*S-CDOT  ) 

X= . 0 1 745*R*PH I 
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00  29  N  = 1  *  NAL 
BN  =  N 

IF  (BN*DEL<1>  -  SOOT)  29«30.30 

29  CONTINUE 
PRINT  7 
GO  TO  200 

30  OF ( N ) =0FL  <  N ) - ( SOOT- <  N- 1 ) *DEL <  1  )  > 

I  NO ( N ) = 1 

I F ( N*EO. 1 ) GO  TO  31 
DE ( N— 1 )=0.0 

31  K*K+ 1 
IFfK.LF.NADGO  TO  14 

TW(K)=TW(K)+OELT*( ATC1K-1 )/<ARHO(K-l )*ACP(K-1 )*DE(K-1 )**2)*(TW{K 
1-1 ) —  T  W (K ) )-ATC(K)/<ARH0<K)*ACP(K)*DEL(K)**2>*<TW(K)-TW(K+l  )  ) ) 
K*K+1 

TW(K)=TW(K)+0ELT*<ATC(K-1 )/(ARHO(K-l )*ACP(K-1 )*DEL(K-1 )**2)* 

1 ( TW ( K— 1 )-TW(< ) ) ) 

PRINT  5  « T I  ME  « ALT , VEL ♦ ANR « T AW  « TEMP , HE . OOOTC , WT . TM , TG  » SDOT , R . 

1 (TWIN) (N=1  « NP2 ) 

«  FORMAT (  AS  NECESSARY  ) 

1F(ALT)200,200, 100 
200  CONTINUE 
END 

0265 
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HEATAB  3,  Version  of  HEATAB 

HEATAB  3,  a  modification  and  extension  of  the  basic  HEATAB  program,  calcu¬ 
lates  the  total  ablation  mass  loss  rate  for  a  conical  reentry  vehicle.  The 
body  is  divided  into  a  number  of  segments,  or  stations,  and  the  heating  rate 
at  each  station  is  computed.  In  addition  to  the  input  data  required  in  HEATAB, 
the  fore  and  aft  radii  and  slant  lengths  of  each  body  segment  must  be  supplied 
to  the  computer.  This  input  into  the  program  is  identified  by  the  comment 
statement  "Additional  Vehicle  Data."  Computation  of  the  skin  temperature 
profile  is  identical  to  the  method  used  in  HEATAB,  except  that  this  calculation 
is  made  for  each  x-station  for  each  trajectory  point. 

To  compute  the  total  mass  loss  rate,  the  following  procedure  is  used.  The 
weight  loss  rate  per  unit  area  is  first  computed  then  the  instantaneous  and 
total  recessions  are  calculated.  Using  these  numbers  and  the  initial  radii 
and  slant  length  of  the  particular  segment  of  interest,  the  mass  loss  rate 
per  uni»-  area  is  multiplied  by  the  appropriate  surface  area  to  give  a  local 
integrated  mass  loss  rate.  Summing  these  local  rates  over  the  body  results  in 
the  total  weight  loss  rate. 

As  written,  the  HEATAB  3  output  will  give,  for  each  time  interval,  an 
output  of  trajectory  conditions  and  body  conditions.  The  trajectory  conditions 
are  items  such  as  altitude,  velocity,  and  Mach  number.  Body  conditions  listed 
for  each  x-station  are  Reynolds  number,  weight  loss  rate  per  unit  area,  wall 
recession,  and  wall  temperature. 
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TWCI.J) 

IND(I,J) 

JIND(I) 

VCP(I,J) 

VTC(I,J) 

V8H0(I,J) 

DEL(I, J) 

DE(I.J) 

SDOT 

XI  (J) 

RS(J) 

RL(J) 

SL(J) 

WT(J) 

TM(J) 

TMELT 

RE 

DELT 

TEMP 

TAW 

SIGMA 

GVIS 

AKL1 

AKL2 

AKT1 

AKT2 

CHA 

NOFX 

TIME 

ALT 

AMACH 

VELS 

RHO 


LIST  OF  SYMBOLS,  HEATAB  3 

Layer  temperature  at  lamina  I  and  station  J 

0  *  virgin  material 
1  *  char  material 

0  •  not  ablating 
1  ■  ablating 

Virgin  specific  heat 

Virgin  thermal  conductivity 

Virgin  density 

Original  lamina  thickness 

Instantaneous  lamina  thickness 

Total  recession  of  wall 

Location  of  x-station,  ft 

Smaller  radius  of  cone  segment  for  a  particular  x-station 
Larger  radius  of  cone  segment  for  a  particular  x-station 
Slant  height  at  station  X 

Instantaneous  weight  loss  rate  at  station  X,  lb/sec 
Total  weight  flux  at  station  X,  lb/ft2-sec 
Ablation  temperature,  8R 
Transition  Reynolds  number 
Time  step,  sec 

Ambient  temperature  at  altitude,  ALT 
Adiabatic  wall  temperature,  9R 
Stephan-Boltrmann  constant 
Ablation  gas  viscosity 
(See  KEATAB  Symbol  Listing) 

(See  HEATAB  Symbol  Listing) 

Cone  half-angle 
Number  of  x-stations 
Time,  sec 
Altitude,  ft 
Mach  number 

Acoustic  velocity,  ft/sec 
Ambient  density,  slugs/ft3 
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VIS 

SPC 

VEL 

Q 

TCA 

VISG 

PE 

CPA 

PBL 

P 

GCP 

GTC 

PR 

ANR 

RF 

HE 

DTR 

TG 

GE 

ALPG 

QDOTC 

TWK 

VT 

VC 

CC 

CT 

QDOTR 

DMD 

PXE 

TMDOT 

CDOT 

WDOT 


Ambient  viscosity 

Cone  surface  pressure  coefficient 

Velocity,  ft/sec 

Dynamic  pressure  lb/ft2 

Air  thermal  conductivity 

Air  viscosity 

Ambient  pressure,  lb/ft2 

Specific  heat  of  air 

Boundary  layer  edge  pressure,  lb/ft2 

Boundary  layer  pressure,  atms 

Ablation  gas  specific  heat 

Ablation  gas  thermal  conductivity 

Prandtl  number 

Reynolds  number 

Recovery  factor 

Convective  heat  transfer  coefficient 

Dynamic  temperature  rise 

Gas  temperature 

Gas  emissivity 

Absorbtivity  at  gas 

Convective  heat  flux 

Dummy  variable  used  in  subroutines 

Dummy  variable  (virgin  thermal  conductivity) 

Dummy  variable  (virgin  specific  heat) 

Dummy  variable  (char  specific  heat) 

Dummy  variable  (char  thermal  conductivity) 

Hot  gas  radiation  heat  llux 

Diffusion  controlled  oxxdscion  mass  loss 

Dummy  variable  in  TFiDOT  equation 

Mass  loss  rates  in  sublimation  region 

instantaneous  recession 

Integrated  weight  loss  rate,  lb/sec 
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HEATAB  IXI 


JOYOOOl  ,6, 100  .50000, YOUNG. 5791  .WLDE.2241  , 

WBK,  75000,  ,,,,30000* 

HFATA03, 


PQ OGRA*-  HFATAB3( INPUT. OUTPUT) 

AEPOOYNAMJC  HEATING  FOP  R/v 

THIS  MODIFICATION  OF  HEATAB  GIVES  THF  ABLATION  MASS  LOSS  RATE 
FOR  A  RE-ENTERING  BOOY.  ^ 

DIMENSION  TW<6« 1 0) « I  NO (6. 10 ) «  JINDI 10 ) t VCP  <6. 10) «VTC(6. 1 0 ) «VRHO<  6. 
110) .DEL (6. 10) «0E(6« 10) .SDOT ( 1 0 ) . X 1 < 1 0 ) .RS < 1 0 ) »RL {  1 0 ) . SL ( 1 0 ) . WT { 1 0 ) 
2.TM(10> 

CONSTANTS  DEFINED 
Pt-3. 14  150 
P=2.71B 
TMELT  *6760. 

RE  *  l.E+7 
DELTa.05 
TEMP=2Rfl.2 
TAW=TEMP 


G  =  32.160 
AJ  =  77R. 

SIC-MA  =  4.8F-13 
GVIS  a  4.42QE-5 
AKL 1  =  5370. 

AKT1  =  4240, 

AKL2  a  5.37 
AKT2  *  «,77 
PRINT  BO  1 
BO  1  FORMAT (1  HI) 

0  IDENTIFICATION  CARD 

READ  502 
50?  PORMATIROH 
1 

PRINT  502 

r  INPUT  -  CONE  HALF  ANGLE.  NUMBER  OF  X-STATIONS.  X-STATIONS 

READ  503.  CHA.NOFX* (XI t I ) ♦  Jsl.NOFX) 

50?  FORMAT (F10. 4, By, I  5 ,5F 1 O . 4/BF  10.4) 

PRINT  BOA , CHA , NOFX . (X! (I >.1*1 «NOFX) 

BOA  FORMAT ( 2 <, 1 BHCONE  HALF  ANGLE  =  .F8.2.5X ♦ 23HNUMBER  Oc  X-STATIONS  = 
1  ,  I  5  ,  /  ,  5X  ,  i  OHX«\STAT  I  ONS  ♦  /  .  (5X.F10.3)  ) 

DO  1  J  a  1 . NOFX 
J I ND { J )  =  0 
VCP ( 5 . J  ) ,  a  .315 
VCR(6,J)  a  ,208 
VTC (F, J) x,7AE-04 
VT C ( 6 «  J )  =  .03694 
VRHO ( 5 ♦ J )  a  91,7 
VRHO ( 6 «  J )  a  169. 

00  1  I  a  1,6 
TKI(IiJ)  a  5A4, 

1  I ND ( I . J )  a  0 

C  ADDITIONAL  VEHICLE  DATA 

PFAD  50b ,  ( PL ( J  )  . Ja 1  , NOFX ) 

PF AD  505,  <DS(J) ,Jal  , NOFX) 

READ  505. (SLCJ) «Jaj , NOFX) 

5CE  FORMAT (8F1 0.5) 

DO  o?  I  =  1,3 
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DO  95  J= 1 t NOFX 
DEL(I.J)  =  <  0083? 

R5  r>F  (  I  ,  J)  =  .^ar 33 
A'O  Q6  I  =  1  ,NO^X 
orL(4,n  =  «oco4i 
DP  (4«t)  =  .00041 
DFL (5,1)  =  .00030 
DP  (5.1)  =  .00333 
npL (6.1)  =  .005 
96  DP  (6.1)  =  .005 

".  INPUT  -  TIWF.  ALTITUDE.  AND  MACH  NUMBFP  FPOM  TRAJFCTOPY 

500  PFAD  506. TIME, ALT, AMACH 

506  FORMAT (OF 10.2) 

IF  ( T I Mp-20 .50 )  1000,1000.403 
1000  CALL  VFLCF ( VELF , ALT )  J 

CALL  TF'VPA  (  TEMP  «  ALT  ) 

CALL  RH0F(PH0, ALT) 

CALL  VISCO( VIS, ALT) 

CALL  9>PCV  (  ?PC  ,  AMACH  ) 

VEL=VELS*AMACH 
0  =  .5*RH0*V£L**2 

517  PRINT  507,  TIME  .ALT,  VEL,  AMACH,  0 

507  FORMAT (/.2X.7HTI ME  =  , F7 . 2 , 5X . 1 1  HALT  1 TUDE  =  . F 1 0 • 2 , 5X *  1 1 HVELOC I TY 
1=  .F9.2.5X, 10HMACH  NO  a  .F5.2.5X. 1 9HDYNAM I C  PRESSURE  a  ,E12.6./, 
22X.7HSTATI0N.6X. 15HREYN0LDS  NUMBER , 5X , 26HWF I GHT  LOSS  RATE  (LBS/SEC 
3) ,5X, 15HT0TAL  RECESS ! ON , 5X » 1 6HWALL  TEMPFRATURE) 

CALCULATF  FOP  EACH  X-STATJON 
DO  401  LeN  =  1 , NOFX 

X  =  XI (LFN) 

IF  (SOOT (LFN)-. 00833)  6,7,7 
>  K=  1 


GO  TO  15 
I  F  <  5D0T ( LFN ) 
K  =  2 

I ND ( < , LFN )  = 
GO  TO  15 
I  F ( SOOT ( LFN ) 


-  .01 666 ) B , 9 , 9 


-  .02499)  10,11,11 


OVEPHFAT* ) 


IND(K.LFN)  =  1 
GO  TO  15 

11  IF(SDOTfLFN)  -  .0341)  12,13,13 

K  =  4 

!ND(K,LFN)  =  1 
GO  TO  15 
PP 1  NT  14 

FORMAT <*  OVEPHFAT*) 

GO  TO  4C3 

15  TP  =  TPMP+ , 58* ( Tw ( K , LEN ) -TEMP )  +  , 1 9* ( TAW-TEMP ) 
CALL  TCAT(TCA.TR) 

CALL  V  I  FOOT (VICG,TR) 

PF=21 16.?*F**(-4,25509F-05*ALT) 

CD A  = , 

PBL.  =  «°C*0+PE 
P  =  PFJL/?1  16.2 

I F ( J I ND (LFN) .FO.O) GO  TO  16 
CALL  GCPP(GCP.P) 

GTC=4 .4?OF-05* ( GCP+.297 ) 


dq=GCo*gv!5/GTC 
GO  TO  17 

16  PRsCP. -*V  1 FG/TCA 

17  ANRsVPL*X*OHO/V!S 
I F ( ANR-PE )  1  8  «  1 0  , 1° 

18  RF=PP**,5 
HF.=  .575#TCA/X#PR**  •  33*ANR**  .  5 
GO  TO  20 

IQ  RF=PP**,33 

HF= • 02®6*TC A/X*PR** • 33*ANR+* , 0 

20  DTR*VEL**2./<2,*G*AJ*CPA) 

TAW=TEMP+RF*DTR 
TG=TEMP+DTR 
IFCTG-8500. >22,21 ,21 

21  GE=. 02169 
GO  TO  23 

22  GF=-.085*F*#(-.61  144E-03*TG) 

23  ALPG*  «  OP* { \  G/TMELT )  **«45 
ODOTC  «  HE*(TAw-TW(K,LEN)  ) *DELT 
IF< IND(K,LEN) ,EQ. 1 >GO  TO  26 

24  IF ( TW { K *LFN) -TMELT  >25,26,26 

25  TWK=TW(K ,LEN) 

CALL  VTCPC ( VT , TWK ) 

CALL  VC.OpC  <  VC  ,  TWK  ) 

VTC ( < *LpN )  =VT 
VC.P  <K,LFN>  «  VC 
VRHQ(K,LFN)*90,4 
IF(K.EQ, 1 )G0  TO  27 

TW  <  K , LEN )»TW(K, LEN ) +DELT*  <  VTC ( K- : »LEN ) / ( VRHO ( K- 1  , LEN ) * vCP ( K- I , L 
1*DE<K-1 , LEN )**2)*(TW<K-1 ,LEN ) -TW ( K , LEN ) ) -VTC ( K ,LEN ) /  ( VRHO(K,LE 
2*VCP<K,LFN)*DE(K,LEN)**2)*<TW<K,LEN)-  TW(K+1 ,LEN> ) ) 

50  t  F  (  TW ( K  ,LtrN) -  TMFLT )  36,51  ,51 

51  TW(K,LFN)=TMELT-1 • 

GO  TO  36 

?6  TWK  =  TW(K,LEN> 

CALL  CTCPC(C.T,TWK> 

CALL  CCPP ( CC , TWK ) 

VCP ( K , LEN )  =  CC 
VTC 1 K ,LEN ) =  CT 
VPHO(K,LcN) =74, 

27  TW(K,'.FN)*TW(K,LEN)+DELT/{VCP(K,LEN)*VRHO(K,LEN)*DE(K,LEN)  )  *  (  WE  "f  (  T 
1 AW-T W  <  K , Lc N ) )-VTC(K,LEN)/DE(K,LEN)*(TW(K,LFN>-TW(K+l .LEN)  ) ) 

20  I F ( TW ( K  » LFN ) -TmELT >72 ♦ 29 ,29 

72  J  t  NO ( LPN )  =  0 
GO  TO  36 

29  TW(K,LEN) =TMELT 
I  NO (K ,LFM ) = 1 
J  t  NO  v  LEN )  =  1 
C  ALLGC.R0  {  GCP ,  P  ) 

OOOTRsE 1 GMA* ( GF*TG**4-ALPG*TMELT**4 ) *DELT 
I F ( OOOTP ) 3 1 ,31 ,3? 

31  ODOTRsC, 

32  I F ( ANR— PE  >33,33,34 

33  OMD= ( GDOTC+ODOTR ) / ( AKL 1 + AKL2* ( GCP*TAW-CPA*TMELT ) ) 

GO  TO  35 

34  DMD=<0P0TC  +  QD0TR)/<AET1+AKT?*(GCP*TAW-CPA*TMFLT>  j 

?«  dxc-1  1  ,OesF+04/TMELT 


40 


2  n 


T'-'DCT  =  r'MD*(l.  +  2.  645+0^/  ( PRL** .67*F**PXf  )) 

TM  (  LFM  )  sTM  <  LEN )  +Tf/nn  x 
~ ^OT ( LFN ) =TM ( LFN ) /74 • 

CnOT-TMDOT/74. 

U«7  (LEN)=»I*(PL  (LFM )  +oc  <  LCN  ) •  *SDGT  (  LFN  )  )  *C00T*90 . 4*  SL(LEN) 
JTND(Lcm)=1 

iir  (  c.^ot  (LFN)-.CFP,3?  )  “?*tr'4,  F4 

*  54  I F< SHOT (LFN)-.oi 666)5^,56, 56 

56  Ic'(.r'00T(LFN)-.r>?/»oo)f:7,cB,Kn 

fR  DE< 4,LcN)=DrL{4,LEN)-(SD0T<LEN)-3.*DEL<  1  «  LFN )  ) 

*  GO  TO  36 

57  OF  (  3 , LEN ) =0FL {  3  .LFN  )  -  (  S.DOT  (  LEN )  -  2  •*DEL  (  1  *  LEN )  ) 

00  TO  *»6 

c5  DF ( ? « L~  N ) * DEL <  ?  « LEN ) - ( SDOT { LFN ) -DEL ( 1 ♦ LEN )  ) 

GO  TO  *>6 

05(1 ,L5M)-OCL( 1 ,LEN) -FOOT < LEN > 

^  K  =  K+1 

I F ( K  »  GT  #4 ) GO  TO  37 
GO  TO  34 

37  TW ( 5  • LFN ) =TW ( 5 »LEN ) +DELT * ( VTC ( 4  »LEN ) / <  VRHO ( 4  « LEN ) *VCP ( 4  »LEN ) *DE ( 4  « 
1  LFN  )  **?  )  *  <  TW  <  4  ,LEN  )  -  TW  (  5 ,  LEN )  )  -VTC  (  5  *  LEN )  /  (  VRHO  <  5  ,LEN )  *VCP  (  5  .  LEN  )  * 
ED^L  (  5  » LCN  )  **?  )  * ( TW ( 5  * LFN  )  — T W  <  6  ♦  LEN  1  )  ) 

TW(6,L'rN)  *TW(6,LFN)+DFLT*(  VTC(5»LEN)/{  VRHO  (  5 ,  LEN )  *VC.P  <  5  ,  LEN )  *OFL  (  5 
1 ,L“N)**?)-(TW(P,LEN)-TW(6,LEN) ) ) 

41  POJNT  509 « y . ANP  * WT ( LEN ) «  SDOT ( LEN ) , TW ( 1  * LEN ) 

500  FORMAT  <F) 0,3,5V, F 14,8 ♦ 1 4V , F9 . 4 , 20X . F7. 5 . 1  OX  * F6 . 0 ) 

401  CONTINUE 

WOOT  =  0.0 

C  CALCULATE  TOTAL  WEIGHT  LO^S  RATE 

00  40?  LL  =  1  » NOFX 

*  WOOT  a  WOOT  +  WT(LL) 

40?  CONTINUE 

POINT  FOR, ALT, VFL. WOOT 

SCR  FORMAT (3X.19HFOR  AN  ALTITUDE  OF  « F 1 0 . 2 . 23HFEET  AND  A  VELOCITY  OF  , 
1  Ffl  »  2  » 33HFT/.SIT  C  «  THF  WEIGHT  LOSS  RATF  IS  <  F£  #  3  *  9HLBS/SEC  •  > 

GO  TO  500 
403  CONTI  Ml  IP 
END 
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